Namn och CID på gruppmedlemmar: Blend Ahmed Omar (blend), Albin Östling (ostlinga), Ebbe Ledin (ebbel)

import numpy as np # Standard paket för att hantera matamatik och arrayer
import matplotlib.pyplot as plt # Standard paket för att plotta figurer
plt.style.use("ggplot")
import scipy.io as io # Vi lånar funktionen 'io' från scipy för att smidigt kunna importera .mat-filer

# Funktioner för HUPP:en

def fft2c(x):
    '''
    2D Fourier transform
    
    Denna är perfekt som den är. Bara att använda!
    '''
    return np.fft.fftshift(np.fft.fft2(np.fft.fftshift(x)))

def ifft2c(x):
    '''
    2D inverse Fourier transform
    
    Denna är perfekt som den är. Bara att använda!
    '''
    return np.fft.fftshift(np.fft.ifft2(np.fft.fftshift(x)))

Uppgift 1 - Skriv en Python-funktion som implementerar PAS!

Utgå från skelletkoden nedan. Funktionen är ej klar! Byt ut alla '.x.' för att få PAS funktionen att fungera.

def PAS(E1, L, N, a, lam0, n):
    '''
    Funktion för att propagera E1 sträckan L genom PAS
    '''
    
    # Varje sampelpunkt i k-planet motsvarar en plan våg med en viss riktning [kx,ky,kz]
    delta_k =2*np.pi/(N*a)                                            # Samplingsavstånd i k-planet
    
    kx      = np.arange(-(N/2)*delta_k, (N/2)*delta_k, delta_k) # Vektor med samplingspunkter i kx-led
    ky      = kx                                                # och ky-led
    
    KX, KY  = np.meshgrid(kx,ky)                                # k-vektorns x- resp y-komponent i varje 
                                                                # sampelpunkt i k-planet
    
    k =2*np.pi*n/lam0                                             # k-vektorns längd (skalär) för en plan våg i ett material med brytningsindex n
    
    KZ = np.sqrt(k**2-KX**2-KY**2, dtype=complex)                   # k-vektorns z-komponent i varje sampelpunkt.
                                                       # dtype=complex tillåter np.sqrt att evaluera till ett komplext tal
    
    fasfaktor_propagation = np.exp(1j*KZ*L) # Faktor för varje sampelpunkt i k-planet
                                           # multas med för att propagera sträckan L i z-led 

    A  = a**2/(2*np.pi)**2*fft2c(E1)                # Planvågsspektrum i Plan 1
    B  = A*fasfaktor_propagation        # Planvågsspektrum i Plan 2 (Planvågsspektrum i Plan 1 multat med fasfaktorn för propagation)
    E2 = N**2*delta_k**2*ifft2c(B)
    
    return E2

Uppgift 2 - Gauss bleibt Gauss (gammalt tyskt ordspråk: Gauss förblir Gauss)

Kolla att en gaussisk stråle förblir gaussisk vid propagation (bara $𝜔$ ändras) ända till fjärrfältet genom att kolla hur fältet ser ut på några olika avstånd fram till $𝐿=𝑓$ (lämplig fokallängd för linsen före Plan 1 kan vara $𝑓=10$ cm och $1/e^2$-radien på infallande fält $𝜔_{in}=1$ mm, se PDF för $𝜔$-definition). Behöver ej redovisas!

# Propagera fält med PAS
# GLÖM EJ ATT ÄNDRA VARIABLER MED '.x.'

N               = 2**10                 # NxN är antalet samplade punkter (rekommenderad storlek N=1024)
sidlaengd_Plan1 = 4e-3                  # Det samplade områdets storlek (i x- eller y-led) i Plan 1 (rekommenderad storlek 4 mm)
a               = sidlaengd_Plan1/N     # Samplingsavstånd i Plan 1 (och Plan 2 eftersom vi använder PAS)
L               = 100e-3                # Propagationssträcka (dvs avstånd mellan Plan 1 och 2)

lambda_noll = 633e-9                    # Vakuumvåglängd för rött ljus från en HeNe-laser
n_medium    = 1                         # Brytningsindex för medium mellan Plan 1 och 2
k           = 2*np.pi*n_medium/lambda_noll                     # K-vektorns längd 

# Definera koordianter i plan 1
x = np.arange(-(N/2)*a, (N/2)*a, a)     # Vektor med sampelpositioner i x-led
y = x                                   # och y-led

X, Y = np.meshgrid(x, y)                # Koordinatmatriser med x- och y-värdet i varje sampelposition
R    = np.sqrt(X**2 + Y**2)             # Avståndet till origo för varje sampelpunkt

# Definera lins och cirkulär aperatur
f_lins = 100e-3                         # Fokallängd på linsen före Plan 1
T_lins = np.exp(-1j*k*R**2/(2*f_lins))  # Transmissionsfunktion för en lins (linsen är TOK)

D_aperture = 2e-3                       # Diameter för apertur
T_aperture = R < (D_aperture/2)         # Transmissionsfunktion för en cirkulär apertur ("pupill")

# Definera fält i plan 1
# GLÖM EJ ATT ÄNDRA VARIABLER MED '.x.'

omega1      = 1e-3                      # 1/e2-radie (för intensiteten, dvs 1/e-radie för amplituden) för infallande Gaussiskt fält
E1_in_gauss = np.exp(-R**2/omega1**2)   # Infallande fält: Gaussiskt med plana vågfronter och normalinfall (dvs konstant fas, här=0)
E1_in_konst = np.ones(X.shape)          # Infallande fält: Konstant i hela plan 1. np.ones(X.shape) ger en matris fylld med ettor som har samma storlek som X


E1_gauss    = E1_in_gauss*T_lins         # Fältet i Plan 1 (precis efter linsen) för gaussisk stråle 
E1_cirkular = E1_in_konst* T_lins* T_aperture  # Fältet i Plan 1 (precis efter linsen) för konstant fält som passerat genom cirkulär apertur


E1          = E1_gauss              # Välj fall!

I1      = np.abs(E1)**2     # Intensiteten är prop mot kvadraten på fältets amplitud (normalt struntar man i proportionalitetskonstanten)
I1_norm = I1/np.max(I1)     # Då vi inte är intresserade av det absoluta värdet av intensiteten så kan det vara trevligt att noramlizera intensiten innan vi plottar

# Plotta intensitet i plan 1 
x_mm = x*1e3
y_mm = y*1e3

plt.figure()
image = plt.imshow(I1_norm, extent = [x_mm.min(), x_mm.max(), y_mm.min(), y_mm.max()])
plt.colorbar(image)

plt.title(r'Intensitet i plan 1. Verkar OK, eller?')
plt.xlabel(r'x $[$mm$]$')
plt.ylabel(r'y $[$mm$]$')
plt.grid(alpha=0.2)

# Plotta fas i plan 1
plt.figure()
image = plt.imshow(np.angle(E1), extent = [x_mm.min(), x_mm.max(), y_mm.min(), y_mm.max()])
plt.colorbar(image)

plt.title(r'Fas i plan 1. Verkar OK, eller?')
plt.xlabel(r'x $[$mm$]$')
plt.ylabel(r'y $[$mm$]$')
plt.grid(alpha=0.2)

# Propagera till plan 2
# GLÖM EJ ATT ÄNDRA VARIABLER MED '.x.'
L = 1e-2

E2      = PAS(E1, L, N, a, lambda_noll, n_medium)         # Propagera med vår PAS funktion
I2      = np.abs(E2)**2    # Intesitet i plan 2
I2_norm = I2/np.max(I2)    # Normaliserad intensitet i plan 2 för att plotta
# Plotta intensitet i plan 2 
plt.figure()
image = plt.imshow(I2_norm, extent = [x_mm.min(), x_mm.max(), y_mm.min(), y_mm.max()])
plt.colorbar(image)

plt.title(r'Intensitet i plan 2 för L = 1 cm. Verkar OK, eller?')
plt.xlabel(r'x $[$mm$]$')
plt.ylabel(r'y $[$mm$]$')
plt.grid(alpha=0.2)
# Plotta fas i plan 2
plt.figure()
image = plt.imshow(np.angle(E2), extent = [x_mm.min(), x_mm.max(), y_mm.min(), y_mm.max()])
plt.colorbar(image)

plt.title(r'Fas i plan 2 för L = 1 cm. Verkar OK, eller?')
plt.xlabel(r'x $[$mm$]$')
plt.ylabel(r'y $[$mm$]$')
plt.grid(alpha=0.2)

L = 5e-2

E2      = PAS(E1, L, N, a, lambda_noll, n_medium)         # Propagera med vår PAS funktion
I2      = np.abs(E2)**2    # Intesitet i plan 2
I2_norm = I2/np.max(I2)    # Normaliserad intensitet i plan 2 för att plotta
# Plotta intensitet i plan 2 
plt.figure()
image = plt.imshow(I2_norm, extent = [x_mm.min(), x_mm.max(), y_mm.min(), y_mm.max()])
plt.colorbar(image)

plt.title(r'Intensitet i plan 2 för L = 5 cm. Verkar OK, eller?')
plt.xlabel(r'x $[$mm$]$')
plt.ylabel(r'y $[$mm$]$')
plt.grid(alpha=0.2)
# Plotta fas i plan 2
plt.figure()
image = plt.imshow(np.angle(E2), extent = [x_mm.min(), x_mm.max(), y_mm.min(), y_mm.max()])
plt.colorbar(image)

plt.title(r'Fas i plan 2 för L = 5 cm. Verkar OK, eller?')
plt.xlabel(r'x $[$mm$]$')
plt.ylabel(r'y $[$mm$]$')
plt.grid(alpha=0.2)

L = 8e-2

E2      = PAS(E1, L, N, a, lambda_noll, n_medium)         # Propagera med vår PAS funktion
I2      = np.abs(E2)**2    # Intesitet i plan 2
I2_norm = I2/np.max(I2)    # Normaliserad intensitet i plan 2 för att plotta
# Plotta intensitet i plan 2 
plt.figure()
image = plt.imshow(I2_norm, extent = [x_mm.min(), x_mm.max(), y_mm.min(), y_mm.max()])
plt.colorbar(image)

plt.title(r'Intensitet i plan 2 för L = 8 cm. Verkar OK, eller?')
plt.xlabel(r'x $[$mm$]$')
plt.ylabel(r'y $[$mm$]$')
plt.grid(alpha=0.2)
# Plotta fas i plan 2
plt.figure()
image = plt.imshow(np.angle(E2), extent = [x_mm.min(), x_mm.max(), y_mm.min(), y_mm.max()])
plt.colorbar(image)

plt.title(r'Fas i plan 2 för L = 8 cm. Verkar OK, eller?')
plt.xlabel(r'x $[$mm$]$')
plt.ylabel(r'y $[$mm$]$')
plt.grid(alpha=0.2)

L = 10e-2

E2      = PAS(E1, L, N, a, lambda_noll, n_medium)         # Propagera med vår PAS funktion
I2      = np.abs(E2)**2    # Intesitet i plan 2
I2_norm = I2/np.max(I2)    # Normaliserad intensitet i plan 2 för att plotta
# Plotta intensitet i plan 2 
plt.figure()
image = plt.imshow(I2_norm, extent = [x_mm.min(), x_mm.max(), y_mm.min(), y_mm.max()])
plt.colorbar(image)

plt.title(r'Intensitet i plan 2 för L = 10 cm. Verkar OK, eller?')
plt.xlabel(r'x $[$mm$]$')
plt.ylabel(r'y $[$mm$]$')
plt.grid(alpha=0.2)
# Plotta fas i plan 2
plt.figure()
image = plt.imshow(np.angle(E2), extent = [x_mm.min(), x_mm.max(), y_mm.min(), y_mm.max()])
plt.colorbar(image)

plt.title(r'Fas i plan 2 för L = 10 cm. Verkar OK, eller?')
plt.xlabel(r'x $[$mm$]$')
plt.ylabel(r'y $[$mm$]$')
plt.grid(alpha=0.2)

Uppgift 3 - Gauss-strålens minsta spotsize

Kolla tumregeln om minsta spotsize, $D_{spot}=𝐶 \frac{𝜆}{D_{start}} 𝐿$, där $𝐶≈1$. Gör det för specialfallet gaussisk stråle ($𝜔_{in}=1$ mm) som fokuseras med hjälp av lins med $𝑓=10$ cm respektive $𝑓=1$ m. För definition av stråldiametern $𝐷$, använd i båda planen $𝐷=2𝜔$, där $𝜔$ är $1/e^2$-radien för intensiteten i respektive plan.

# KOD
# Propagera till plan 2
# GLÖM EJ ATT ÄNDRA VARIABLER MED '.x.'

f_1 = 10e-2
f_2 = 1

def D_spot_gauss(f):
    
    f_lins = f                        # Fokallängd på linsen före Plan 1
    T_lins = np.exp(-1j*k*R**2/(2*f_lins))
    
    E1    = E1_in_gauss*T_lins

    E2      = PAS(E1, f, N, a, lambda_noll, n_medium)         # Propagera med vår PAS funktion
    I2      = np.abs(E2)**2    # Intesitet i plan 2
    I2_norm = I2/np.max(I2)
    I2_norm_flat = I2_norm.flatten()  #Vi flattenar matriserna så att vi enkelt kan hitte det R där I2_norm är närmast 1/e^2
    R_flat = R.flatten()
    
    index = np.argmin(np.abs(I2_norm_flat - 1/np.exp(2)))
    omega = R_flat[index]
    return 2*omega

C = D_spot_gauss(f_1) *2e-3 /(633e-9*f_1)  # Vi löser ut C ur tumregelekvationen
print(C)

1.258644232225837

# Plotta intensitet i plan 2 
plt.figure()
image = plt.imshow(I2_norm, extent = [x_mm.min(), x_mm.max(), y_mm.min(), y_mm.max()])
plt.colorbar(image)

plt.title(r'Intensitet i plan 2. Verkar OK, eller?')
plt.xlabel(r'x $[$mm$]$')
plt.ylabel(r'y $[$mm$]$')
plt.grid(alpha=0.2)

C = D_spot_gauss(f_2) *2e-3 /(633e-9 *f_2)
print(C)

1.2835703001579646

f_lins = f_2                        # Fokallängd på linsen före Plan 1
T_lins = np.exp(-1j*k*R**2/(2*f_lins))
    
E1    = E1_in_gauss*T_lins

E2      = PAS(E1, f_2, N, a, lambda_noll, n_medium)         # Propagera med vår PAS funktion
I2      = np.abs(E2)**2    # Intesitet i plan 2
I2_norm = I2/np.max(I2)


plt.figure()
image = plt.imshow(I2_norm, extent = [x_mm.min(), x_mm.max(), y_mm.min(), y_mm.max()])
plt.colorbar(image)
plt.title(r'Intensitet i plan 2. Verkar OK, eller?')
plt.xlabel(r'x $[$mm$]$')
plt.ylabel(r'y $[$mm$]$')
plt.grid(alpha=0.2)

Vad blir det mer exakta värdet på $𝐶$ för den gaussiska strålen med den valda definitionen av $𝐷$?

1.258644232225837 för f = 10 cm och 1.2835703001579646 för f = 1m så vi skulle säga att ett mer exakt värde skulle ligga vid 1.27

Uppgift 4 - Cirkulära strålens minsta spotsize

Gör samma sak som i uppgift 3 fast nu ska fältet i Plan 1 ha konstant intensitet i ett cirkulärt tvärsnitt (diameter $D_{start}$), och $D_{spot}$ i Plan 2 definieras som ”innersta mörka ringens” diameter (vilket är en ganska generös definition av spotdiameter, så $𝐶$ bör bli klart större än 1).

def D_spot_cirk(f):
    
    f_lins = f                        # Fokallängd på linsen före Plan 1
    T_lins = np.exp(-1j*k*R**2/(2*f_lins))
    
    E1    = E1_in_konst*T_lins*T_aperture

    E2      = PAS(E1, f, N, a, lambda_noll, n_medium)         # Propagera med vår PAS funktion
    I2      = np.abs(E2)**2    # Intesitet i plan 2
    I2_norm = I2/np.max(I2)
    I2_norm_flat = I2_norm.flatten()  #Vi flattenar matriserna så att vi enkelt kan hitte det R där I2_norm är närmast 1/e^2
    R_flat = R.flatten()
    
    index = np.argmin(np.abs(I2_norm_flat - 1/np.exp(2)))
    omega = R_flat[index]
    return 2*omega

C1 = D_spot_cirk(f_1) *2e-3 /(lambda_noll*f_1)
print(C1)
C2 = D_spot_cirk(f_2) *2e-3 /(lambda_noll*f_2)
print(C2)

1.6558560259922914
1.6272758416123274

Vad blir det mer exakta värdet på $𝐶$ för den cirkulära strålen och de valda definitionerna av $𝐷$?

1.6558560259922914 för f = 10cm och 1.6272758416123274 för f = 1m

f_lins = f_2                        # Fokallängd på linsen före Plan 1
T_lins = np.exp(-1j*k*R**2/(2*f_lins))
    
E1    = E1_in_konst*T_lins*T_aperture

E2      = PAS(E1, f_2, N, a, lambda_noll, n_medium)         # Propagera med vår PAS funktion
I2      = np.abs(E2)**2    # Intesitet i plan 2
I2_norm = I2/np.max(I2)


plt.plot(Y[:, 0], I2_norm[:, I2_norm.shape[1]//2])  # Intensitetsprofil längs en linje
plt.xlabel("Y-position (m)")
plt.ylabel("Intensitet")
plt.title("Intensitet i fjärrfältet för f = 1m")
plt.grid(True)
plt.show()

f_lins = f_1                        # Fokallängd på linsen före Plan 1
T_lins = np.exp(-1j*k*R**2/(2*f_lins))
    
E1    = E1_in_konst*T_lins*T_aperture

E2      = PAS(E1, f_1, N, a, lambda_noll, n_medium)         # Propagera med vår PAS funktion
I2      = np.abs(E2)**2    # Intesitet i plan 2
I2_norm = I2/np.max(I2)

plt.xlim(-0.0002, 0.0002)
plt.plot(Y[:, 0], I2_norm[:, I2_norm.shape[1]//2])  # Intensitetsprofil längs en linje
plt.xlabel("Y-position (m)")
plt.ylabel("Intensitet")
plt.title("Intensitet i fjärrfältet för f = 10 cm")
plt.grid(True)
plt.show()

Kan du i något av fallen $𝑓= 10$ cm eller $𝑓=1$ m se svaga tecken på numeriska fel i simuleringen av intensitetsfördelningen i fjärrfältet?

graferna visar att intensitetsfördelningen är mycket brusigare för f = 1m och vi tror att detta kan vara ett tecken. Detta skulle även vara logiskt för alla fält som "åker" långa sträckor i ett homogent fält kan expandera vilket betyder att de kan expandera tillräckligt för att punkter ska hamna utanför det numeriska fönstret, alltså N b x N b.

Uppgift 5 - Möt en annan medlem i gauss-strålens storfamilj!

Det finns en hel lycklig familj med fält som delar gausstrålens egenskap att inte ändra form (förutom storlek) när de propagerar: så kallade hermite-gaussiska strålar. En av dessa strålar fås helt enkelt genom att multiplicera den (vanliga) gaussiska strålen i startplanet med $𝑥$. Kolla att detta infält också håller formen vid propagation ända till fjärrfältet! Behöver ej redovisas!

T_lins2 = np.exp(-1j*k*R**2/(2*f_2)) 

E1_hermitgauss = E1_in_gauss*X

E2      = PAS(E1_hermitgauss, 10e-2, N, a, lambda_noll, n_medium)         # Propagera med vår PAS funktion
I2      = np.abs(E2)**2    # Intesitet i plan 2
I2_norm = I2/np.max(I2)  

plt.figure()
image = plt.imshow(I2_norm, extent = [x_mm.min(), x_mm.max(), y_mm.min(), y_mm.max()])
plt.colorbar(image)

plt.title(r'Intensitet i plan 2. Verkar OK, eller?')
plt.xlabel(r'x $[$mm$]$')
plt.ylabel(r'y $[$mm$]$')
plt.grid(alpha=0.2)

plt.figure()
image = plt.imshow(np.angle(E2), extent = [x_mm.min(), x_mm.max(), y_mm.min(), y_mm.max()])
plt.colorbar(image)

plt.title(r'Fas i plan 2. Verkar OK, eller?')
plt.xlabel(r'x $[$mm$]$')
plt.ylabel(r'y $[$mm$]$')
plt.grid(alpha=0.2)

#Vi har kollat med många fokallängder mellan 0 och 1 meter och formen bibehålls.

Uppgift 6 - Dubbla budskap

# Ladda in DOE såhär!
DOE = np.flip(io.loadmat('T_DOE_gen2.mat'))["T_DOE_gen2"]

print(DOE)

[[0 0 0 ... 0 0 0]
 [0 0 0 ... 0 0 0]
 [0 0 0 ... 0 0 0]
 ...
 [0 0 0 ... 0 0 0]
 [0 0 0 ... 0 0 0]
 [0 0 0 ... 0 0 0]]

# Skapa det elektriska fältet som träffar ögat och propagera med pas
T_lins3 = np.exp(-1j*k*R**2/(2*0.02))
E1_in = E1_in_konst* T_lins3*DOE
E2 = PAS(E1_in, 0.02, N, a, 633e-9, 1)

# Fältet som har propagerats genom DOE och 'de vises lins' är smidigt att
# plotta i logaritmisk skala för att bilden ska vara tydlig.
# Alltså plotta fältet i plan 2 så här!
I2      = np.abs(E2)**2    
I2_norm = np.log(I2/np.max(I2))  # Log av den normaliserade intensiteten i plan 2

x_mm = x*1e3
y_mm = y*1e3

plt.figure()
image = plt.imshow(I2_norm, extent = [x_mm.min(), x_mm.max(), y_mm.min(), y_mm.max()])
plt.colorbar(image)

<matplotlib.colorbar.Colorbar at 0x159c892e0d0>

T_lins_vise = np.exp(-1j*k*R**2/(2*0.15)) # Vi testade en del styrkor mellan 0 och 0.5, 0.15 var bäst!
E1_in = E1_in_konst* T_lins3*DOE*T_lins_vise
E2 = PAS(E1_in, 0.02, N, a, 633e-9, 1)

I2      = np.abs(E2)**2    
I2_norm = np.log(I2/np.max(I2))  # Log av den normaliserade intensiteten i plan 2

x_mm = x*1e3
y_mm = y*1e3

plt.figure()
image = plt.imshow(I2_norm, extent = [x_mm.min(), x_mm.max(), y_mm.min(), y_mm.max()])
plt.colorbar(image)

<matplotlib.colorbar.Colorbar at 0x159c8294750>

Vilket är det ofarliga meddelandet i den bifogade DOE-koden av generation 2-typ?

Make Newton great again!

Och vilken styrka krävs på ”de vises lins” för att läsa det farliga budskapet i samma kod? Testa dig fram med olika värden på styrkan hos de vises lins tills du kan urskilja budskapet.

0.15 var den bästa styrkan och budskapet var <3 Fresnel